library(ggpubr)
library(pheatmap)
library(RColorBrewer)
library(greekLetters)
library(tidyverse)
library(patchwork)
source('DE_cells/scripts/0.DE_filtering_funcs.R')

# read data -------
chen2021 <- read_csv('DE_cells/results/chen2021-crc_de_meta.csv.gz')
ye24 <- read_csv('DE_cells/results/AH_cOME_de_meta.csv.gz')

# DE cell frac for tumor and normal tissue -----
chen2021_dec <- chen2021 |>
  calc_de_cell_frac(tissue) |>
  mutate(type = case_when(str_detect(type, 'DE') ~ type,
                          .default = str_glue('{type}-')) |>
           fct_relevel("Ahmed DE", "Full DE", 'CD79A-', 'CD79B-', 'CD79-', 'CD247-', 'CD3D-', 'CD3E-', 'CD3G-', 'CD3-')) |> 
  arrange(type) 

crc_chen_pval <- chen2021_dec |>
  select(-fraction) |>
  filter(tissue %in% c('Normal', 'Primary')) |>
  mutate(total = total - count) |>
  pivot_wider(names_from = tissue, values_from = c(total, count)) |>
  rowwise() |>
  mutate(pval = t.test(c(rep(0,total_Normal), rep(1,count_Normal)),
                       c(rep(0,total_Primary), rep(1,count_Primary)))$p.value,
         psignif = case_when(pval < .001 ~ '***',
                             pval < .01 ~ '**',
                             pval < .05 ~ '*',
                             .default = ''))

chen2021_dec |>
  filter(tissue %in% c('Normal', 'Primary')) |>
  mutate(tissue = case_match(tissue, 'Primary' ~ 'Tumor', .default = 'Normal')) |>
  ggplot(aes(type, log10(fraction+1e-5), color = tissue, group = tissue)) + 
  geom_path() +
  geom_text(data = crc_chen_pval, aes(label = psignif, x = type, y = -1.5), inherit.aes = FALSE) +
  theme_pubr() +
  labs_pubr() +
  expand_limits(y = c(-5, -1)) +
  scale_color_manual(values = c('blue', 'red')) +
  labs(x = 'DE cell defining criteria',
       y = 'log10(DE cell fraction + 1e-5)') +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = .5))

crc_cd79 <- chen2021_dec |>
  select(-fraction) |>
  filter(tissue %in% c('Normal', 'Primary')) |>
  mutate(total = total - count) |>
  pivot_wider(names_from = tissue, values_from = c(total, count)) |>
  rowwise() |>
  mutate(pval = t.test(c(rep(0,total_Normal), rep(1,count_Normal)),
                       c(rep(0,total_Primary), rep(1,count_Primary)))$p.value ,
         normal_es = t.test(c(rep(0,total_Normal), rep(1,count_Normal))) |> broom::tidy(),
         tumor_es = t.test(c(rep(0,total_Primary), rep(1,count_Primary))) |> broom::tidy()) |>
  filter(type == 'Full DE')

normal_es <- crc_cd79 |> pull(normal_es) |>
  select(estimate, conf.high)

tumor_es <- crc_cd79 |> pull(tumor_es) |>
  select(estimate, conf.high)

bind_rows(list(normal = normal_es, tumor = tumor_es), .id = 'tissue') |>
  ggplot(aes(tissue, estimate, color = tissue, ymin = estimate, ymax = conf.high)) +
  geom_col(fill = 'white', linewidth = 2) +
  geom_errorbar(width = .3) +
  theme_pubr() +
  labs(y = 'DE cell fraction',
       title = 'DE cell fraction in CRC tumor and normal tissue') +
  scale_color_manual(values = c('blue','red'))

chen2021_dec |>
  write_csv('DE_cells/results/chen2021crc_de_cell_tumor-normal.csv')

## CRC DE per sample -----
crc_de_sample <- chen2021 |>
  filter(tissue %in% c('Normal', 'Primary')) |>
  unite(sub_sample,c(sample, tissue), sep = '.') |>
  calc_de_cell_frac(sub_sample) |>
  separate(sub_sample, into = c('sample', 'tissue'), sep = '\\.') |>
  mutate(type = case_when(str_detect(type, 'DE') ~ type,
                          .default = str_glue('{type}-')) |>
           fct_relevel("Ahmed DE", "Full DE", 'CD79A-', 'CD79B-', 'CD79-', 'CD247-', 'CD3D-', 'CD3E-', 'CD3G-', 'CD3-')) |> 
  arrange(type) 

crc_de_sample |>
  ggplot(aes(tissue, fraction, color = tissue)) +
  stat_summary(fun = 'mean', geom = 'col', fill = 'white') +
  geom_jitter(height = 0, width = .1) +
  stat_compare_means(comparisons = list(c('Normal',)),method = 't.test') +
  scale_y_continuous(expand = expansion(mult = c(0,.3))) +
  scale_color_manual(values = c('red','blue')) +
  theme_pubr() +
  facet_wrap(~type, scales = 'free')

# AH-cOME ------
## cOME & ctrl -----
ye24_dec <- ye24 |>
  calc_de_cell_frac(cOME) |>
  mutate(type = case_when(str_detect(type, 'DE') ~ type,
                          .default = str_glue('{type}-')) |>
           fct_relevel("Full DE", "Ahmed DE", 'CD79A-', 'CD79B-', 'CD79-', 'CD247-', 'CD3D-', 'CD3E-', 'CD3G-', 'CD3-')) |> 
  arrange(type) 

ye24.come_pval <- ye24_dec |>
  select(-fraction) |>
  mutate(total = total - count) |>
  pivot_wider(names_from = cOME, values_from = c(total, count)) |>
  rowwise() |>
  mutate(pval = t.test(c(rep(0,total_ctrl), rep(1,count_ctrl)),
                       c(rep(0,total_cOME), rep(1,count_cOME)))$p.value,
         psignif = case_when(pval < .001 ~ '***',
                             pval < .01 ~ '**',
                             pval < .05 ~ '*',
                             .default = ''))
ye24_dec |>
  filter(type != 'CD3-') |>
  ggplot(aes(cOME, fraction, fill = cOME)) +
  geom_col() +
  facet_wrap(~type, scales = 'free_y') +
  theme_pubr() +
  labs(x = 'Group', title = 'DE cell fraction in AH patient adenoids') +
  stat_pvalue_manual(ye24.custpval, inherit.aes = F) +
  scale_y_continuous(expand = expansion(mult = c(NA, .2)))

ye24.custpval <- ye24.come_pval |>
  filter(type != 'CD3-') |>
  mutate(fraction = max(count_cOME / total_cOME, count_ctrl/ total_ctrl)) |>
  customize_pvalue(y = fraction, p.value = pval, facets = type,
                   group.1 = 'cOME', group.2 = 'ctrl')

ye24_dec |>
  ggplot(aes(type, log10(fraction+1e-5), color = cOME, group = cOME)) + 
  geom_path() +
  geom_text(data = ye24.come_pval, aes(label = psignif, x = type, y = -1.5), inherit.aes = FALSE) +
  theme_pubr() +
  labs_pubr() +
  #expand_limits(y = c(-5, -1)) +
  scale_color_manual(values = c('blue', 'red')) +
  labs(x = 'DE cell defining criteria',
       y = 'log10(DE cell fraction + 1e-5)') +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = .5))

### per sample -----
ah.meta <- ye24 |>
  select(orig.ident, sex, `FCGR2B-I232T`, `IGHG1-G396R`, cOME, pid) |>
  distinct(orig.ident, .keep_all = T)

ah.meta

ye24_per.sample <- ye24 |>
  calc_de_cell_frac(orig.ident) |>
  mutate(id = orig.ident) |>
  left_join(ah.meta)

ye24_per.sample |>
  filter(type != 'Ahmed DE') |>
  mutate(cOME = fct_relevel(cOME, 'ctrl'),
         type = ifelse(type == 'Full DE', 'Full DE', type) |>
           fct_relevel('Full DE')) |>
  ggplot(aes(cOME, fraction, color = cOME)) +
  stat_summary(fun = 'mean', geom = 'col', fill = 'white') +
  geom_jitter(height = 0, width = .1) +
  stat_compare_means(method = 't.test', color = 'black') +
  scale_y_continuous(expand = expansion(mult = c(0,.3))) +
  scale_color_manual(values = c('blue','red')) +
  theme_pubr() +
  facet_wrap(~type, scales = 'free')

## G396R ------
ye24.g2r <- ye24 |>
  calc_de_cell_frac(`IGHG1-G396R`) |>
  mutate(type = case_when(str_detect(type, 'DE') ~ type,
                          .default = str_glue('{type}-')) |>
           fct_relevel("Full DE", "Ahmed DE", 'CD79A-', 'CD79B-', 'CD79-', 'CD247-', 'CD3D-', 'CD3E-', 'CD3G-', 'CD3-')) |> 
  arrange(type)

ye24.g2r.pval <- ye24.g2r |>
  select(-fraction) |>
  mutate(total = total - count) |>
  pivot_wider(names_from = `IGHG1-G396R`, values_from = c(total, count)) |>
  rowwise() |>
  mutate(pval = t.test(c(rep(0,total_GG), rep(1,count_GG)),
                       c(rep(0,total_RR), rep(1,count_RR)))$p.value,
         psignif = case_when(pval < .001 ~ '***',
                             pval < .01 ~ '**',
                             pval < .05 ~ '*',
                             .default = ''))
  
ye24.g2r |>
  filter(type != 'CD3-') |>
  ggplot(aes(`IGHG1-G396R`, fraction, fill = `IGHG1-G396R`)) +
  geom_col() +
  facet_wrap(~type, scales = 'free_y') +
  theme_pubr()

ye24.g2r |>
  ggplot(aes(type, log10(fraction+1e-5),
             color = `IGHG1-G396R`, group = `IGHG1-G396R`)) + 
  geom_path() +
  theme_pubr(x.text.angle = 45) +
  labs_pubr() +
  expand_limits(y = c(-5, -1)) +
  labs(x = 'DE cell defining criteria',
       y = 'log10(DE cell fraction + 1e-5)')

# I232T ------
ye24.i2t <- ye24 |>
  calc_de_cell_frac(`FCGR2B-I232T`) |>
  mutate(type = case_when(str_detect(type, 'DE') ~ type,
                          .default = str_glue('{type}-')) |>
           fct_relevel("Full DE", "Ahmed DE", 'CD79A-', 'CD79B-', 'CD79-', 'CD247-', 'CD3D-', 'CD3E-', 'CD3G-', 'CD3-')) |> 
  arrange(type)

ye24.i2t.pval <- ye24.i2t |>
  select(-fraction) |>
  mutate(total = total - count) |>
  pivot_wider(names_from = `FCGR2B-I232T`, values_from = c(total, count)) |>
  rowwise() |>
  mutate(pval = t.test(c(rep(0,total_II), rep(1,count_II)),
                       c(rep(0,total_IT), rep(1,count_IT)))$p.value,
         psignif = case_when(pval < .001 ~ '***',
                             pval < .01 ~ '**',
                             pval < .05 ~ '*',
                             .default = ''))

ye24.i2t |>
  filter(type != 'CD3-') |>
  ggplot(aes(`FCGR2B-I232T`, fraction, fill = `FCGR2B-I232T`)) +
  geom_col() +
  facet_wrap(~type, scales = 'free_y') +
  theme_pubr()

ye24.g2r |>
  ggplot(aes(type, log10(fraction+1e-5),
             color = `IGHG1-G396R`, group = `IGHG1-G396R`)) + 
  geom_path() +
  theme_pubr(x.text.angle = 45) +
  labs_pubr() +
  expand_limits(y = c(-5, -1)) +
  labs(x = 'DE cell defining criteria',
       y = 'log10(DE cell fraction + 1e-5)')
